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1 .  INTRODUCTION 


Bus  report  discusses  the  numerical  evaluation  of  certain  types  of 
multiple  integrals  that  arise  in  the  theoretical  description  of  short- 
range  coherent  laser  sensing  systems.  This  work  was  part  of  an  effort 
at  the  Harry  Diamond  Laboratories  (HDL)  to  explore  the  fuzing  potential 
of  such  sensors,  which  have  promise  for  significantly  increased 
detection  sensitivity  compared  with  conventional  direct  detection 
systems  and  would  make  Doppler  information  available  to  the  fuze 
designer. 

A  critical  problem  for  the  performance  evaluation  of  coherent  laser 
fuzes  can  be  seen  from  the  following  basic  considerations.  When  a 
nearly  monochromatic  laser  beam  illuminates  a  nominally  flat  target 
surface  at  normal  incidence,  the  reflected  surfaces  of  constant  phase 
immediately  in  front  of  the  target  will  be  nominally  planar  and  repli¬ 
cate  the  surface  roughness  features,  to  the  extent  that  shadowing 
effects  do  not  complicate  the  picture.  If  the  roughness  features  have  a 
scale  comparable  to  the  illuminating  laser's  wavelength — a  situation 
which  holds  for  aircraft  skins  and  the  CO2  laser  wavelength — then  the 
variation  of  phase  of  the  reflected  field  in  a  plane  immediately  in 
front  of  the  target  would  be  expected  to  oscillate  with  amplitude  on  the 
order  of  2n  in  a  random-looking  manner,  as  the  observation  point 
traverses  the  illuminated  region.  The  reflected  field  retains  the 
temporal  coherence  of  the  transmitter  laser,  but  the  spatial  coherence 
is  severely  degraded.  Fortunately,  as  the  reflected  wavefront 
propagates  in  space  toward  a  receiver  system,  spatial  coherence  tends  to 
improve;  that  is,  the  spatial  coherence  within  a  given  receiver  aperture 
increases  with  the  distance  between  the  receiver  and  the  reflector. 
This  improvement  is  important  because  coherent  detection  depends  on 
establishing  a  phase  match  in  space  between  the  reflected  and  local 
oscillator  wavefronts  in  the  mixing  region  at  the  receiver. 


In  coherent  optical  radar  applications,  the  receiver  apertures  and 
target-to-receiver  distances  are  such  that  there  is  virtually  no  spatial 
coherence  degradation  effect.  For  fuzing  applications,  however,  target- 
to-receiver  distances  are  usually  relatively  small--on  the  order  of  10 
m — so  that,  for  practical  receiver  apertures,  spatial  coherence  degra¬ 
dation  effects  must  be  considered.  This  report  describes  a  computa¬ 
tional  apparatus  for  analyzing  these  effects  and  presents  sample 
calculations.  Many  aspects  of  the  computational  problem  present  severe 
numerical  difficulties.  Methods  of  overcoming  these  difficulties  are 
discussed,  and  a  listing  of  the  FORTRAN  source  code  which  carries  out 
the  calculations  is  presented  in  appendix  A. 

Hie  problem  considered  is  the  numerical  calculation  of  the  mutual 
intensity  function  of  the  field  reflected  from  a  rough  surface  at 
arbitrary  distances  from  the  reflector.  Specifically,  we  are  interested 
in  the  mutual  intensity  function  Jb(p^  ,P2)  between  two  points, 
P^  and  P^,  on  an  observation  plane — the  plane  of  the  receiver 
aperture— given  the  mutual  intensity  function  J(J(p1  ,P2)  between  arbi¬ 
trary  points  (P|  ,£2)  i-n  the  illuminated  target  region,  a.  The  mutual 
intensity  function,  defined  precisely  below,  has  the  information  needed 
to  determine  spatial  coherence  degradation  effects  in  the  receiver,  and 
its  propagation  from  the  target  region  can  be  formulated  in  terms  of 
Huygens-Fresnel  integral  propagators.1  In  addition,  the  effect  of  the 
rough  surface  can  be  modeled  through  the  specification  of  the  mutual 
intensity  function  boundary  values  in  the  illuminated  target  region. 


Hie  mutual  intensity  function  is  used  to  describe  the  coherence  and 
intensity  pattern  of  light  in  space.  In  this  sense  it  is  a  generali¬ 
sation  of  the  familiar  concept  of  intensity.  Formally,  the  mutual 
intensity  is  defined  between  two  points,  P1  and  P2,  as  the  long-time 
average  of  the  product  of  the  complex  analytic  signal  of  the  field 
1#T.  Born  and  E,  Wolf,  Principles  of  Optics,  Ch  10,  Pergamon  Press 


t 


at  P1  with  the  complex  conjugate  of  the  analytic  signal  at  £2*  where  the 
analytic  signals  are  evaluated  at  the  same  time.  A  further  generali¬ 
zation  would  be  to  introduce  a  time  element  and  define  a  function  being 
the  time  average  product  of  the  analytic  signals  at  P1  and  P2  when 
these  signals  are  observed  a  set  time  apart.  This  function  is  the  well- 
known  mutual  coherence  function  of  physical  optics.  (Neither  of  these 
functions  preserves  any  polarization  information,  which  would  be  a 
further  generalization.  In  this  report,  however,  we  will  be  exclusively 
concerned  with  the  mutual  intensity  function. ) 

We  assume  that  the  mutual  intensity  function  is  known  on  some  plane 
in  space  (the  target  plane)  and  that  this  plane  is  parallel  to  the 
observation  plane.  In  addition,  we  assume  that  the  mutual  intensity 
function  on  the  target  plane  is  nonzero  only  over  a  circular  area,  a, 
and  within  a  the  function  depends  only  on  the  distance  between  the 
evaluation  points.  That  is,  given  any  two  points,  P1  and  P2,  on  the 
target  plane,  the  mutual  intensity  between  P1  and  P2  is  zero  if  either 
P1  or  P2  is  outside  o;  otherwise,  Ja(Pi»l?2)  =  J0( !?i  "  PjO*  These 
assumptions,  while  somewhat  restrictive,  still  include  the  important 
case  of  a  homogeneously  rough  surface  with  no  preferred  direction 
illuminated  by  a  circular  beam  of  coherent  light. 


2.  PROBLEM  FORMULATION 

We  begin  with  the  general  equation  for  the  mutual  intensity  on  an 
observation  plane  B.  Figure  1  shows  the  geometry  being  used.  We  have1 


ik(si-s2) 


S1S2 


d?1d?2 


(1) 


1M.  Born  and  E.  Wolf,  Principles  of  Optics,  Ch  10,  Pergamon  Press 
(1970). 
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Figure  1.  Geometry  of  mutual  intensity  propagation  problem.  (Line  00' 
is  normal  to  both  a  and  B. ) 

In  the  above,  k  is  the  mean  wavenumber,*  Sj  =  |Pj  -  Pj | ,  i  =  /-I ,  and  Aj 
are  directionality  factors  given  by 

Aj  =  t  cos  h  ' 

where  is  the  angle  between  (Pj  -  Pj )  and  the  normal  to  0,  and  \  is 
the  wavelength.  (The  integration  is  limited  to  0  because  we  assume 
Jo(?1,?2)  zero  outside  0.)  From  figure  1  we  have 

■  (2) 

*The  field  is  assumed  to  be  nearly  monochromatic . 


so  we  obtain 


JB(Si'?S)  -  (t?)2  Jo  Jo 


ik(si-s2) 


d!ld?2 


(3) 


The  equations  in  this  form  could,  in  theory,  be  integrated  numeri¬ 
cally  simply  by  defining  a  coordinate  system  and  performing  the  fourfold 
nested  integration  indicated.  Typically,  however,  the  wavenumber  k  i 
quite  large  and,  unless  one  is  willing  to  restrict  attention  to  lorn 
range  calculations,  the  quantity  k(s1  -  s2)  in  the  exponential  create 
severe  efficiency  problems.  Additionally,  J0(pifp2)  typically 
finely  structured  function,  with  a  scale  of  about  1 0-l+  mm.  For  even 
moderate-size  target  o,  this  implies  excessively  large  numbers  of  grid 
points  for  straightforward  integration. 

The  above  formulation  can  be  simplified  within  the  assumption  pre¬ 
sented  above.  Since  Jc(p.,,P2)  depends  only  on  |p1  -  P2|,  we  can 
approximately  decompose  it  into  a  constant  term  (also  termed  the 
residual  contribution)  and  a  sum  of  terms  corresponding  to  varying-size 
circles  of  constant-value  mutual  intensity.  That  is,  we  can  determine 
appropriate  Wj,  a..,  and  n  such  that,  for  a  given  error  e. 


n 

Z 

3=o 


w3Ja 


(4) 


where 


a 

J0(x,a) 


x  <  a 

V 


(5) 


The  residual  term  is  obtained  by  setting  aQ  =  <».  (We  note  that 
is  9enerally  complex,  so  in  general  the  Wj  are  complex.)  We 
thus  need  to  consider  only  the  integration 


5b(?1'£2)  =  ®)2  la  la  “  ?2l'a) 


ik(si-s2) 

"•^2 .  d?1d?2 


for  an  arbitrary  a.  This  integration  will  be  considered  in  the  next 
section. 


3.  CALCULATION  OF  INTEGRAL 


We  are  now  concerned  with  the  integral 


VEi'Ei)  -  ($?)*  fo  fo  -  ?2 1 ' rc ) 


ik(s1-s2) 


dP.dP. 


(where  we  have  dropped  the  circumflex  notation  for  and  J0  used  in  the 
last  section  and  used  rc,  rather  than  a,  as  more  expressive  of  a  "cutoff 
radius").  The  behavior  of  J  (|Pj  -  p2l'rc)  can  *3e  accounted  for  by 
changing  the  limits  of  integration 


JB(?i'?2)  =  (5S) 2  la  X(p2)  —37^~  d?id?2  ' 


where  o(p2)  is  the  intersection  of  a  with  a  circle  of  radius  rc  centered 
on  P2»  This  integral,  in  turn,  splits  naturally  into  two  terms: 


( 1 )  a  "central  contribution"  term 


where  o  -  o'  is  the  central  portion  of  o  with  radius  r_  -  r_  (where  r_ 
is  the  radius  of  or)  and 

(2)  an  "edge  effect"  term. 


ik(Sl-s2) 


a2s2 
S1  s2 


dP1dP2 


(10) 


where  o'  is  the  remaining  annulus  of  breadth  rc  on  the  outside  edge  of 

o.  The  purpose  of  this  partition  of  the  integral  is  to  isolate  the 

first  term,  where  o(p2]  does  not  intersect  the  edge  of  a  (and  is  thus  a 

circle),  from  the  second  term,  where  a(P2)  is  guaranteed  to  intersect 

the  edge  of  o,  and  hence  is  a  more  complicated  case.  As  we  shall  see 

later,  realistic  problems  correspond  to  r  <<  r_  (except  for  the  resid- 

u  o 

ual  term),  and  hence  a'  is  a  relatively  small  portion  of  a. 


We  note  that  if  r  >  r  the  central  contribution  does  not  exist; 
c  -  a 

thus,  the  only  contribution  is  due  to  the  "edge  effect,"  where  a’  =  ff. 
Additionally,  if  rc  >  2ra,  the  intersection  of  a(P2)  with  a  is  a  itself 
(since  P2  is  inside  a),  so  we  obtain  in  this  case 


jb(Ei-E2)  -  J2(Ei'?2)  -  ($S)2  So  So 


ik(s1-s2) 


S2S2 
S1  S2 


d?1d?2' 


rc  i 


2r 


(11) 
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l* 


This  is  the  formula  for  the  residual  term  contribution  mentioned  in  the 
last  section.  The  integral  in  this  case  can  be  further  reduced  to  the 
product  of  two  essentially  identical  integrals  as  follows. 


We  have 


=  ®)2 


L  U 


lk(s,-s2) 


TIP~  d?1d?2'  rc  i  2t0 

S1S2 


-  rkDi2 


(t  ■?-)(/• 


(12) 


-iks- 


dP- 


°  s?  -2 


We  now  define  a  function. 


t 


iks 


j(P'l2)  -  Z  d?'  3  =  I?  "  S’ I 


(13) 


for  arbitrary  observation  point  p’  and  arbitrary  circle  Z.  Figure  2 
shows  the  geometry  relevant  to  the  definition  of  j(p'  |z).  We  also 
define  the  notational  convention  that  a  circumflex  over  a  "prime"  point 
denotes  the  projection  of  that  point  onto  the  target  plane.  That  is,  in 
a  coordinate  system  with  the  origin  at  0  in  figure  2,  if 

P'  =  (x',y’,D)  , 


then 


P'  *  (x',y’,0) 


J 

rJ 


3 
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Figure  2.  Geometry  for  calculating  J(p'|e) 


In  terms  of  the  j(p* (e)  function,  we  have 


jb(?v£2)  -  (®2  rc  *  2ro  '  (14) 

where  J*  indicates  the  complex  conjugate  of  J. 

We  now  consider  the  numerical  evaluation  of  J(P’|o).  We  first  note 
that  the  integrand  is  constant  for  all  points  on  any  given  circle  cen- 

A 

tered  at  Pf.  This  is  because  s  is  the  slant  height  of  a  right  circular 

cone  with  apex  at  P’  and  a  base  with  a  circumference  passing  through  P. 

Thus,  the  most  natural  way  to  integrate  over  o  is  to  use  a  polar  coordi- 

* 

nate  system  centered  at  P'.  We  show  this  coordinate  system  in  figure  3, 
where  we  have  arbitrarily  selected  the  0=0  line  (using  an  r  -  0  polar 
coordinate  notation)  to  pass  through  the  center  of  0.  In  this  coordi¬ 
nate  system,  we  have  s  «  +  r^,  and  the  integral  becomes 


13 


j(r ! 


"  -k 


ik/D2+r2 


|0olrJ 

J-eo(r)  d9 


+s  -ik/r^+r* 

0  2r0  (r )  - -  dr 

°  D2  +  r2 


where  0Q(r)  is  the  angle  at  which  the  circle  of  radius  r  intersects  the 


boundary  of  or,  and  r 


S  -  rc,  if  P'  is  outside  o  (fig.  3(a)),  or  rQ  = 


0  if  P'  is  inside  a  (fig.  3(b)).  The  value  of  2r60(r)  is  simply  the  arc 

a 

length  of  the  portion  of  the  circle  centered  at  P'  of  radius  r  falling 
inside  a.  From  simple  geometry,  assuming  that  the  circles  intersect  as 
in  figure  3(a),  we  have 


0Q  =  cos 


A2  * s2  - 

v  J5? — /  ■ 


In  figure  3(b),  if  r  <  r^  -  S,  the  circles  do  not  intersect,  and  we  have 
0Q  =  it;  otherwise,  0Q  is  still  given  as  above. 


We  thus  have  two  cases.  If  P*  falls  outside  a,  we  have 


,  ,  ,  fS+ra  eik/52^2 

J(P*  a)  =  2  /  r  - 

JS-r_  D2  +  r 


cos 


-1 


r2  +  S2  -  r2 


"3Sf" 


-^dr,  S  >  ra  . 


If  P'  falls  inside  a,  we  have 


f  ,  feo-S  eik^ 

j(*‘|o)  =  2lT  Jo  riTTT 


D2+r2 


dr 


(17) 


f^+ra  ^ik/D2+r2  / r2  +•  S2  -  r2\ 

+  2  J  rg-S  r  ~D2~T72~  COB  \ - 5Sr - ‘J  ^ 


(18) 


S  <  r0  . 


The  first  term  of  the  second  case  (S  <  rg)  can  be  integrated  in 
closed  form  in  terms  of  exponential  integrals.  We  are  considering  the 
integral 
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2ir 


r 


eiklA>2+r2 

D2  +  r2 


dr 


(19) 


If  we  change  variables  u 


\/b2  +  r2, 


we  obtain 


The  quantity^!  +  [ (r-  -  s)/d]2  will  occur  again,  so  we  define 


■*W. 


(21 ) 


We  also  define 


(22) 


which  will  also  be  used  later.  We  then  have,  from  the  definition2  of 
the  exponential  integral  E1 ( z ) , 


=  2ir[E1(-ikD)  -  E^-ikD^)] 


(23) 


In  terms  of  cosine  and  sine  integrals,  Ci(x)  and  Si(x),  we  have,  from 
the  identity 


E.j  { -ix)  =  -Ci  (x)  -  iSi(x)  +  i 


(24) 


the  expression 


I1  =  2ir[ci (kDi(>)  -  Ci(kD)]  +  i2ir[si(kDi|;)  -  Si(kD)]  .  (25) 

ZM.  Abramowitz  and  I.  4.  Sfcegrun,  Handbook  of  Mathematical  Functions 
with  Formulas,  Graphs,  and  Mathematical  Tables,  NBS  Applied  Math  Series 
55  (1970). 


16 


While  this  expression  could  be  evaluated  using  any  of  various  methods 
for  calculating  the  Si  and  Ci  functions,  we  note  that  normally  kD  is 
very  large  (on  the  order  of  106),  so  special  precautions  must  be  taken. 
One  method  well  suited  to  evaluating  these  functions  for  large  arguments 
is  to  express  Si(x)  and  Ci(x)  in  terms  of  "auxiliary"  functions  f(x)  and 
g(x)  as2 


Si (x)  = -f  -  f(x)  cos(x)  -  g(x)  sin(x)  , 
Ci(x)  =  f(x)  sin(x)  -  g(x)  cos(x)  , 


(26) 


and  use  asymptotic  expansions  for  f(x)  and  g(x)  for  large  arguments. 


f(x,  -l/,  .  2L  +  4L.5L+  .  .  \ 

X  \  X2  x**  x6  ) 

_^1.  +  5L-21+  .  .  .) 

X2  \  X2  X**  X6  / 


g(x) 


(27) 


For  realistic  values  of  kD  and  kD<»,  these  expansions  converge  extremely 
rapidly;  only  three  or  four  terms  are  needed  for  most  calculations.  At 
long  ranges,  however,  when  D  »  r  -  S,  accuracy  may  be  lost  in  computing 
the  differences  involved  in  evaluating  I1 .  If,  in  the  expressions  for 
Si(x)  and  Ci(x)  in  terms  of  f(x),  g(x),  sin(x),  and  cos(x),  we  express  i|» 
as  1  +  ( ip  -  1)  in  the  arguments  for  sin(x)  and  cos(x),  we  can  arrive  at 
the  expression 


2weikD(a  +  iB) 


(28) 


where 


a  =  g(kD)  +  f (kDi|i)  sin  [kD(ip  -  1)]  -  g(kD^)  cos  [kD(<i  -  1)] 


(29) 


0  -  f (kD)  -  f  (kD\p)  cos  [kD(\p  -  1)]  -  g(kDf)  sin  [lcD(ip  -  1)]  . 


2M.  Abramowitz  and  I.  A.  Stegun,  Handbook  of  Mathematical  Functions 
with  Formulas,  Graphs,  and  Mathematical  Tables,  NBS  Applied  Math  Series 
55  (1970). 


The  factor  ^  -  1  in  the  trigonometric  functions  can  be  computed  at  long 
ranges  using  the  approximation 


✓1  +  x  «  1 


(30) 


to  yield 


(31) 


At  short  ranges,  of  course,  the  expression 


♦  “  1 


(32) 


should  be  used. 


We  now  turn  attention  to  the  second  integral  involved  in  j(p'|o) 
(eq  (17)  and  (18)).  This  integral,  in  both  cases,  can  be  written 


(33) 


This  integral  must  be  evaluated  numerically,  since  no  closed-form  solu¬ 
tion  is  known.  This  integral,  and  several  others  to  be  encountered 
later  in  this  report,  will  be  reduced  to  a  canonical  form  f(x)eixdx 
in  order  to  make  use  of  an  efficient  numerical  integration  program 


developed  at  HDL^  .for  integrals  of  this  form*  The  above  integral  can  be 
reduced  to  canonical  form  with  the  substitution  u  =  kV^2  +  r2,  which 
yields 


1 2  =  2 


fkDu  iu  _  P  *  <kD>2  +  (ks>2  ~  (kra)2l 
I  -rr—  COS-1| -  - -  dU  . 

■'kD*  L  2ksVu2  _  (kD)2  J 


At  long  ranges,  the  above  form  for  I2  is  unsuited  for  numerical 
solution.  (Because  w  -  4*  <<  the  integration  interval  is  very  small 
compared  to  the  midpoint  value  of  the  interval,  and  severe  truncation 
errors  will  occur  when  taking  differences. )  A  solution  is  to  develop  a 
long-range  approximate  form  for  I2.  This  can  be  done  by  expanding  the 
radicals  in  equations  (21)  and  (22),  making  a  suitable  change  of  vari¬ 
ables;  that  is. 


VDi|i  =  kD\1  + 


'r  -  S’ 
a 


•  •  •] 


kD  +  Mro  -  s)2  1 


D  »  |ro  -  S|  , 


3 Theodore  H.  Hopp,  A  Routine  for  Numerical  Evaluation  of  Integrals 
with  Oscillation  Behavior ,  Proceedings,  1979  Army  Numerical  Analysis  and 
Computers  Conference,  El  Paso,  TX  ( 1979 )• 


«kD+^(ra  +  s)2h  ^  "  T 


(36 


D  »  ra  +  S  . 

Changing  variables  v  =  u  -  kD,  we  obtain 


I2  -  2e 


ikD 


f  e  coa"1 

v2  +  2kDv  +  (kS)2  -  (krfl)2 

L  v  +  kD  cos 

L 

2kS/  v2  +  2kDv  -1 

dv. 


(37) 


D  »  rQ  +  S  , 


where 


and 


(*g  +  5)2 


(38) 


To  recapitulate,  we  now  have  expressed  the  function  j(p'|o), 
arbitrary  P'  and  circle  o,  as  the  sum  of  two  integrals,  I,  and  l~. 


for 


can  be  evaluated  using  asymptotic  expansions  of  the  exponential  integral 
E.,(z),  and  I2  has  been  expressed  in  canonical  form  for  numerical  evalu¬ 
ation  for  both  the  general  case  and  the  long-range  case. 

Having  established  how  to  evaluate  j(p' |o),  we  can  now  return  to  the 
original  problem  of  calculating  JB(p^,P^).  We  have  (e<l  (8)  to  (10)) 


+  J2(?i'?2)  '  (39) 


which  reduces  to 


for  rc  >  2rQ.  So  far  we  have  not  considered  how  to  evaluate  J2(p1’ ,P^) 
when  rc  <  2rQ  or  how  to  evaluate  J^Pj  ,P^]  at  all.  Evaluating 
J2(p{  (Pj)  for  rc  <  2rg  will  not  be  considered  in  this  report.  As  men¬ 
tioned  before,  for  realistic  problems,  unless  the  residual  contribution 
is  being  considered  (in  which  case  rc  £  2rff),  we  always  have  rc  <<  rg. 
Since  in  this  case  the  area  of  o’  is  very  small  compared  to  the  area  of 
o,  J2(p} *£2) >  being  an  integral  over  o’,  makes  a  negligible  contribution 
to  JB(pj,Pj).  Before  proceeding  with  the  investigation  of  J1(p.j,PJ!),  we 
will  illustrate  by  way  of  an  example  the  calculations  we  are  now  able  to 
perform. 

If  we  consider  a  target,  0,  which  is  perfectly  smooth  (i.e.,  an 
optically  flat  mirror),  illuminated  uniformly  by  coherent  monochromatic 


light,  the  mutual  intensity  at  the  target  is  everywhere  a  real  constant. 
For  simplicity,  we  will  assume  this  constant  to  be  unity.  The  mutual 
intensity  at  an  observation  plane  parallel  to  this  mirror  is  precisely 
jb(  Pj.Pj)  with  rc  =  «»  (i.e.,  the  residual  contribution).  We  note  that 
when  P^  =  P£,  we  have 

jb(?v?i)  =  ©2  lJ(£{l°)l2  *  (41) 

which,  from  the  definition  of  mutual  intensity,  is  the  (real)  intensity 
at  Pj .  Further,  JB(p.j  ,P^)  can  be  interpreted  as  the  normalized  diffrac¬ 
tion  pattern  from  a  circular  aperture  o.  Sample  calculations  were 
performed  for  j(p||o)  for  D  =  10  m,  tg  =  0.05,  A  *  10.6  pm,  and  for  S 
ranging  from  0  to  0.1  m.  This  choice  of  parameters  places  P{  well  into 
the  near-field  region  of  a.  Plots  of  j(p»|oJ  and  of  JB(pJ,Pl)  are 
shown  in  figures  4  and  5.  (The  code  used  to  generate  these  plots  was 
checked  by  moving  Pj  to  the  far-field  region  of  a  to  verify  the  expected 
Fraunhoeffer  diffraction  pattern. ) 


OFFSET  (cm) 


-TAR8ET  SHADOW 


OFFSET  (cm) 

Figure  5.  Intensity  versus  offset  from  target  axis. 

We  now  return  to  the  problem  of  calculating  J.(pj,P£).  Recall  that 
is  given  by  (eq  (9)) 


VEi-ES)  -  (§)2  X(P2)  e  (S^3|Z)  ir-iar-2  ■ 


which  can  be  written 


-  -iks~  T  iks- 

Ji(Ei>Ej)  -  ©  [ /„( p2)  dEi J  "2  • 


But  we  note  that  the  function  j(p'|o)  can  be  introduced,  since  o(p2)  is 


a  circle.  We  have 


-  /0(p2) 


■  (w)2  /o-0-  j[e;W?2)1 


From  its  definition,  j[p.j|a(P2)]  has  circular  symmetry  about  P2«  Equi- 

A 

valently,  it  has  circular  symmetry  about  P^ .  Since  the  major  work  in 
evaluating  the  integrand  of  J^Pj,]?^)  evaluat;*-ng  J[p}  |  o(P2)  ] »  it 

A 

is  natural  to  express  the  integral  in  polar  coordinates  about  P^ .  The 
geometry  involved  is  shown  in  figure  6.  We  define  the  6  =  0  axis  as 

A  A 

passing  through  P}  and  P£.  Defining  the  aperture  function 


we  obtain 


u(r , 0;o)  =  | 


1  if  (r,0)  is  inside  o 
0  otherwise 


r(5i'?2)  “  ©2/ro+(r°"rc)  rJ[?iK?2)] 


[/?  u{r'9; 


-iks  "] 

a  -  a')  — -  d0  dr  , 


I 


where  R1 ,  from  figure  6,  is  the  distance  from  Pj  to  the  center  of  a 
and  rQ  is  given  by 

1" 

0  if  Pj  is  inside  o  -  o' 

R1  -  (r  -  rcJ  if  Pj  is  outside  a  -  o' 


Figure  6.  Geometry  for  calculating  central  contribution 


When  r  <  (r  -  rc)  -  ,  u(r,0;o  -  o')  *  1  for  all  values  of  0. 

A 

This  occurs  only  when  Pj  is  inside  a  -  o'  (as  shown  in  fig.  6).  When 

(r0  -  *c)  _  Ri<  r  <  Ri  +  (ra  “  rc)»  there  is  a  pair  of  values  (0Q,  0.J 

such  that  u(r,0;o  -  o')  =  1  only  when  0Q  <  6  £  0,  (also  shown  in  fig. 

6).  From  the  diagram  we  see  that  there  is  a  0*  such  that  0Q  =  0R  -  0' 

and  01  =  0R  +  0'.  If  Pj  is  at  (xj,y^,o)  and  P£  is  at  (x£,y£,o),  with 
the  origin  at  the  center  of  o»  we  have  from  simple  geometry 

*  -V*i'2  +  yj2  , 

r2  =  Vx$2  +  y'22  , 


(50) 

d  =V(*i  "  x^)2  +  (yj  -  y$J2  , 


.  /Rf  +  d2  -  R§\ 

«R  "  COS  1  \—-"3Rld - l)  ’ 

and 

/r*  +  R?  -  (r  -  t  )A 
6'  -  cos  1  \— - - tlj  . 


When  Pj  falls  outside  o  -  o'  (R.  >  r  -  r  )  we  have 


iK.+ir.-r 

Ji(Ei-Ej)  -  (§)  J,  fr  ,  ,  "[p)|o(p2)J 

J  R1“lro"rcJ 

f  feR+e*  -iks2  1 

I  5 — o —  d0  dr 

L  Jv9'  sl  J 


When  falls  inside  o  -  o'  (R,  <  rQ  -  rc),  we  have 

-  (®2  fn  °  1  rJ[?il®(!2)] 


ftir  "iks2 


d0  dr 


„  rfr-r  1+R, 


l~  i~0R+6'  -iks2 

~ — 5 —  dQ  dr  . 

L  J  0R-0'  32  J 


rJ  [?ild(?2)] 


At  this  point  we  have  two  integrals  to  evaluate  to  obtain  J^(pj,p£)  and 
thence  JB( ,P^),  namely,  the  outer  integral  of  equations  (51)  and  (52), 
which  is  an  integration  in  r,  and  the  inner  integral  in  0,  which  has  the 


two  forms 


From  the  equations  (50)  for  0R  and  0’,  we  have 


0  <  0R,0'  i  *  > 

thus, 

I 

-IT  <  0R  -  0'  <  TT  , 

|  and 

0  <  0R  +  0'  <  2u  . 

I 

We  identify  four  cases: 

(1)  0R  -  0'  >  0,  0R  +  0*  i  it  . 

! 

I 

In  this  case  the  variable  transformation  in  question  has  a  single-valued 
inverse  on  the  range  of  integration. 

(2)  0R  -  0'  <  0,  0R  +  0'  <  rr  . 


In  this  case  the  integration  can  be  treated  as  follows: 


The  last  step  is  taken  to  avoid  integrating  over  [0,8']  twice 


(3)  0R  -  0'  >  0,  0R  +  0'  >  v  . 


The  integration  is  treated  as  follows: 


(using  the  substitution  0  «•  2tt  -  0) 


f*  f' 

I  ...  d0  +  f 
j0D-0'  J  2v-(  0„+6 


30 


•  •  d9  f 

»’) 


V6' 


!L. 


2ir-(6  +6') 


Hie  last  step  follows  from  observing  that  since  0R  <  ir,  0R  -  0'  <  2ir 

(0R  +  0-). 

(4)  0R  -  0*  <  0,  0R  +  0*  >  ir  . 

The  integration  is  treated  as  follows: 

/0R+0-  ro  fir  feR+e  • 

•  •  •  d0  =  I  •  •  •  <30  +  I  •  •  •  <30  +  I  •  •  •  <30 

0P-01  Je  p-0*  Jo  J  7r 


r0,-0R  p  p 

—  I  •  •  •  <30  +  /  •  •  •  <30  +  I 

J  0  JO  J  2tt- 


(0R+0*) 


P'-e  R  f2*- 

=  2  I  ...  d0  +  I 
JO  J0'- 


>(0R+0') 


•  •  •  <30 


f 

J2ir- 


(V01) 


•  •  •  <30  • 


In  the  last  step,  since  0*  <  tt,  we  are  assured  that  0*  -  0R  <  2tt  - 
(0R  +  9').  In  all  cases,  the  integral  in  0  can  be  broken  down  into  a 
sum  of  integrals,  each  of  which  is  over  an  interval  [0L,  0„],  where 


we  have 


<30 


_ 2udu _ 

V(u2  -  A2)(b2  -  u2) 


(65 


u  |  0_e  =Va2  +  2rdk2(l  -  cos  9  ) 


(66 


u|0=0  =  Vb2  -  2rdk2(l  +  cos  0^) 

Finally,  we  obtain 


iks- 


fi 

d0  =  2k2  I 


Vb2 - 2rdk2 [ 1 +cos0y) 


Va2  +2rdk2 ( 1 -COS0L) 


,iu 


du 


iiV(  u2  -  A2  ]  (  b2  -  u2 ) 


(67 


(68) 


As  in  the  integration  for  J2(pj  'Vi)'  when  D  »  r  +  d  numerical  problems 
arise  because  the  interval  of  integration  is  small  compared  to  the 


location  of  the  interval.  In  the  integral,  we  have  two  differences  of 
large  numbers  to  consider,  namely,  u  -  A  and  B  -  u.  If  we  perform  the 
change  of  variable  x  «-  u  -  A,  we  obtain 


pOy  iks 
I  e 


-  d6  =  2k2e“  rVB2-2rdk2(n-cos  0„)  -  A 
^Vft2+2rdk2(  1  -cos  0L)  -  A 


(x  +  A)-yx(B  -  A  -  x )  ( x  +  2A)  ( x  +  A  +  B) 

Now  all  the  differences  between  large  numbers  can  be  calculated  before 
hand.  We  are  concerned  with  the  quantities 


A2  +  2rdk2 ( 1  -  c 


Xy  =  y  B2  -  2rdk2(l  +  cos  0y)  -  l 


B  -  A  . 


Substituting  for  A  and  B,  and  using  the  expansion 


krd(l  +  cos  0y) 


(r  +  d)2  -  rd(l  +  cos  0y) 


\  "2  - 


X 


1 


(r  +  d)2  -  rd( 1 
D2 


(79) 


(r  +  d)2[(r  +  d)2  -  2rd(l  +  cos  0y) ] 

8D4 

This  second  substitution,  in  turn,  works  well  unless  the  interval  of 
integration  is  near  the  singularity  at  A.  This  happens  only  when 
|cos  0y  -  cos  0LI  <<  1  and  |1  -  cos  0jJ  <<  1.  Since  the  interval  of 
integration  cannot  be  entirely  located  near  A  and  entirely  located  near 
B  at  the  same  time,  one  or  the  other  of  the  above  substitutions  can  be 
used  for  any  problem.  We  will  use  the  first  substitution,  x  «•  u  -  A, 
unless  cos  0^  <  0.  In  that  event  we  will  use  the  substitution  y  «-  u  - 
B,  since  |1  -  cos  0LI  <<  1  cannot  hold. 

At  this  point  we  have  specified  how  to  do  all  the  calculations 
for  except  the  integral  in  r.  The  integrands  in  this  case 
cannot  be  expressed  in  the  canonical  form  eixf(x)  in  any  useful  way,  so 
another  method  of  integration  is  needed.  As  will  be  shown  in  the  next 
section,  the  behavior  of  the  integrand  in  this  integral  is  quite  varied 
in  character,  depending  strongly  on  the  separation  of  the  observation 
points  P^j  and  P^.  Thus,  the  number  of  points  at  which  v  will  be  neces¬ 
sary  to  evaluate  the  integrand  cannot  be  predicted  wex  Rather  than 
using  a  fixed  Gaussian  quadrature  scheme  where  all  the  work  done  at  one 


3' 


level  of  approximation  is  useless  at  the  next  level  of  accuracy,  it 
would  be  more  sensible  to  use  an  integration  scheme  that  makes  use  of 
all  previous  work  when  trying  to  increase  the  precision  of  the  integra¬ 
tion.  One  of  the  best  such  schemes  is  an  adaptive  Gaussian  scheme 
modified  to  force  all  points  previously  used  to  be  included  at  the  next 
level  of  precision.4 5  Since  such  a  routine  has  already  been  developed, ^ 
we  will  not  elaborate  on  it  here.  Rather,  we  will  proceed  to  describe 
the  results  of  calculations  for  a  sample  problem. 


4.  SAMPLE  PROBLEM — IDEALIZED  GOODMAN  SURFACE 


In  the  previous  section  we  described  how  to  evaluate  the  mutual 
intensity  between  two  points  on  an  observation  plane  when  the  mutual 
intensity  of  the  reflected  radiation  at  the  target  plane  is  given  by 


1  if  Ip,  -  P2|  <  rc  and 

P,  and  P2  are  inside  a  (80) 


0  otherwise 


for  an  arbitrary  cutoff  radius  rp.  The  justification  for  considering 

such  a  simple  model  for  J  was  briefly  mentioned  in  section  1 .  Any  i 

0 

that  depends  only  on  the  distance  between  the  evaluation  points 

P,  and  P2  can  be  approximated  arbitrarily  well  by  a  finite  sum  of 

functions  of  the  above  form  with  various  r  ’s;  moreover,  such  a  J 

c  O 

dependence  is  expected  when  a  flat,  homogeneously  and  isotropically 


4T.  N .  L.  Patterson,  The  Optimum  Addition  of  Points  to  Quadrature 
Formulae,  Math,  Comp.,  22_  (1968),  847-856. 

5T.  N.  L.  Patterson,  Algor ithm  for  Automatic  Numerical  Integration 
Over  a  Finite  Interval,  Algorithm  468  Comm,  of  the  ACM,  16,  11  (November 
1973),  694-699. 
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rough  target  is  illuminated  by  a  coherent  light  beam.  Further  justifi¬ 
cation  for  considering  the  form  of  equation  (80)  is  furnished  by  Good¬ 
man's  recent  discussion6  of  a  simple  physical  model  of  for  flat  rough 
surfaces  whose  slope  distribution  is  concentrated  mainly  in  a  low  slope 
region. 


An  example  of  Goodman's  predictions  for  JQ  is  shown  by  the  solid 
curve  in  figure  i ,  which  assumes  uniform  reflected  intensity  of  unity  at 
a;  also  shown,  by  dotted  lines,  is  a  rough  approximation  to  the  pre¬ 
dicted  Jo~curve  by  two  functions  of  the  form  in  equation  (80).  It  is 
clear  from  figure  7  that  very  good  approximations  of  Goodman-surface 
J  's  can  be  obtained  by  using  only  a  few  different  values  of  r  . 

In  this  section,  we  apply  the  results  of  section  2  to  an  example, 
namely,  that  of  a  target  whose  JQ  corresponds  to  the  Goodman  surface  of 
figure  7.  We  approximate  by  the  dotted  step  function  in  figure  7  and 
calculate  JB(Pj '?2^  at  an  observation  plane  for  various  Pj ,££  pairs. 

1  _ -APPROXIMATED  MUTUAL  INTENSITY  (SEE  TEXT) 


i  i  i  i  I 


CORRELATION  LENOTH 


SEPARATION  (mm) 

Figure  7.  Mutual  intensity  at  Goodman  surface. 

6J.  W.  Goodman,  Statistical  Properties  of  Laser  Speckle  Patterns , 
Laser  Speckle  and  Related  Phenomena,  Ed.,  J.  C.  Dainty,  Topics  in 
Applied  Physics,  9,  Springer-Verlag  (1975),  63-68 . 


In  its  simplest  form,  Goodman's  model  describes  JQ  in  terms  of  the 
second-order  statistics  of  the  surface  height  variations.  Specifically, 
we  have 


=Vi(p i)i(p2)  exp(-°e  \}  ~  PhU?i  -  ?2 1 )])  •  (81) 


In  the  above,  is  the  Goodman  model  J  and  I(P)  is  the  re- 

fleeted  intensity  at  point  P.  In  the  exponential,  a|  is  the  surface 
height  variance  scaled  to  the  effective  wavelength,  given  by 


o|  =  [k(1  +  cos  $)]2o2h  , 


(82) 


where  k  is  the  wave  number,  0  is  the  angle  of  illumination  relative  to 
the  surface  normal,  and  is  the  actual  variance  of  the  surface  height. 
Finally  Ph(x)  is  the  normalized  surface  height  autocorrelation  function, 
which,  for  simplicity,  is  assumed  to  be  of  the  form 


phU?i  -  ?2D  =  exp\- 


I?1  - 


(83) 


where  rQ  is  the  surface  roughness  correlation  length. 


The  distance  r 


is  the  separation  at  which  the  surface  height  correlation  reaches  a 
value  of  1/e.  Figure  7  plots  a  typical  Jq(P-|»P2)  versus  lpi  ~  P2 1  • 


versus 


The  numerical  parameter  values  used  to  produce  figure  7,  and  which 
will  be  used  for  the  sample  problem  of  this  section,  are  as  follows: 

Wavelength  (A)  =  10.6  um, 

Surface  height  standard  deviation  =  2  ym,  (84) 

Angle  of  illumination  (0)  =0, 

Surface  correlation  length  (rQ)  =  1.06  mm  (=  100X). 

These  values  correspond  to  a  surface  that  is  optically  flat  over  dis¬ 
tances  on  the  scale  of  100  wavelengths,  illuminated  perpendicularly  by  a 
10.6-ym  laser  (e.g.,  a  CO2  laser),  and  observed  at  a  distance  of  10  m. 
Hills  and  valleys  on  the  surface  average  2  pm  of  deviation  from  the  mean 
surface  height.  These  values  yield,  for  the  other  quantities  used  in 
the  Goodman  model. 


Og  =  5.6217 


(85) 


ph(x)  =  exp(-889,996  x2)  . 


For  our  sample  problem  we  take 


where  rc  is  chosen  as  that  point  separation  for  which  =  ^1  +  e  ®J/2; 
i.e.,  the  separation  at  which  Jr  has  fallen  halfway  from  its  maximum 
value  of  unity  to  its  minimum  value  e  The  desired  value  of  rc  is 

readily  found  as 


which  for  the  chosen  numerical  values  of  the  sample  problem  gives  rc  = 
3.83445  x  10-4  m.  Equation  (86)  gives  the  dotted  line  approximation  to 
JG  shown  in  figure  7.  To  complete  the  specification  of  our  sample 
problem,  we  take 


Range  (D)  =  10  m 

Target  spot  radius  (rQ)  =  0.05  m  . 


(88) 


(Note  that  these  values  are  the  same  as  were  used  for  the  sample  dif¬ 
fraction  pattern  calculations  of  the  last  section.) 


The  contribution  to  J0(p},p£)  from  the  first  term  of  equation  (86) 


(compare  with  eq  (4)  and  (40)).  We  call  this  part  of  JB  the  residual 
contribution  since  it  arises  from  the  low  dotted  plateau  in  figure  7, 
which  represents  the  residual  phase  correlation  in  the  reflected  field 
at  o  in  the  limit  of  large  point  separation.  The  residual  phase  cor¬ 
relation  is  physically  related  to  the  specular  component  of  reflection 
from  the  target.  Figure  4  plots  the  modulus  and  phase  of  J(P'|o)  appro¬ 
priate  to  our  sample  problem.  The  contribution  of  the  second  term  in 
equation  (86)  (which  corresponds  to  the  high  phase  correlation  shown  by 
JG  for  small  point  separations)  to  JB(p},P^)  is  more  difficult  to  repre¬ 
sent  since  it  depends  not  only  on  the  separation  of  P^  and  P£,  but  also 

a 

on  their  location  within  the  shadow  of  o.  (We  will  denote  by  a  the 
shadow  of  o  on  the  observation  plane.)  We  note  that  JB(p.j,P£)  is  invar¬ 
iant  to  rotation  about  the  center  of  a  since  the  problem  is  circularly 
symmetric.  We  also  note  that  JB(p',p^)  is  invariant  to  reflection 
through  the  center  of  o.  Thus,  one  convenient  way  of  examining 
JB(P',P^)  is  to  fix  one  of  the  points  (say,  P')  at  some  distance  from 
the  center  of  a  and  move  the  other  point  (p^)  away  from  Pj  at  various 
angles  from  the  radius.  The  geometry  of  this  layout  is  shown  in  figure 
8. 


Figures  9  to  11  show  the  results  of  calculations  for  the  sample 
problem  for  three  positions  of  the  fixed  point,  P{.  Figure  9  shows  what 
happens  when  Pj  is  at  the  center  of  o  (R1  =  0)  and  P£  is  moved  away  from 
Pj.  (Since  this  case  is  circularly  symmetric  about  Pj ,  there  is  no 
need  to  plot  what  happens  at  different  angles  from  a  radius.)  As  can  be 
seen,  the  sharp  step  function  character  of  the  approximated  mutual 
intensity  at  the  target  plane  (fig.  7)  has  been  spread  out  sub¬ 
stantially,  although  it  is  still  apparent  that  the  mutual  intensity 
between  Pj  and  P£  is  stronger  when  the  two  points  are  close  together. 


43 


MV* 


ANGLE  FROM  RADIUS  -  I 


Figures  10  and  11  show  the  results  when  R1  =  1  cm  and  =  3  cm, 
respectively.  For  each  of  these  cases,  P^  was  moved  away  from  P j  at 
five  angles,  as  indicated,  from  the  radius  passing  through  Pj.  As  can 
be  seen  from  the  figures,  the  mutual  intensity  magnitude  as  a  function 
of  d  (fig.  8)  does  not  vary  significantly  for  the  various  values  of  0, 
nor  for  the  three  positions  of  Pj .  The  mutual  intensity  phase,  however, 
does  show  significant  variation  from  case  to  case.  In  all  cases,  the 
phase  exhibits  very  sharp  jumps  as  d  increases.  However,  the  variation 
of  the  phase  as  a  function  of  9  is  noticeably  more  pronounced  for  R1  =  3 
cm  (fig.  11)  than  for  R1  =  1  cm  (fig.  10).  Of  course,  the  phase  does 
not  vary  at  all  as  a  function  of  angle  for  R1  =  0,  since  there  is  cir¬ 
cular  symmetry. 

The  characteristics  of  the  mutual  intensity  can  be  described  in 
general  terms  as  follows.  To  a  first  approximation,  the  magnitude  of 
the  mutual  intensity  between  two  points  within  o  varies  with  d  in  the 
manner  shown  in  figure  9,  regardless  of  where  the  two  points  lie  within 
a.  The  phase  of  the  mutual  intensity  varies  with  d  qualitatively,  like 
the  case  when  Pj  is  at  the  center  of  a*  However,  the  phase  is  modulated 
by  some  perturbation  function  of  8  for  nonzero  values  of  R1 .  This 
perturbation  becomes  more  strongly  a  function  of  0  as  P^  is  placed 
farther  out.  The  specific  nature  of  this  perturbation  is  difficult  to 
ascertain,  however,  from  the  calculations  shown  here. 

The  real  part  of  JB  is  of  the  most  direct  interest  for  optical 
heterodyne  receiver  analysis,  which  is  the  goal  of  the  computational 
apparatus  being  discussed  in  this  report.  The  major  components  of  an 
optical  heterodyne  receiver  are  a  photodetector/narrow-band  amplifier 
combination  and  a  local  oscillator.  The  local  oscillator  is  a  highly 
coherent  optical  source  that  illuminates  the  photodetector  at  relatively 
high  power  levels,  whose  frequency  is  slightly  shifted  relative  to  that 


of  the  received  (target  reflected)  radiation.  McGuire7  has  discussed 
the  dependence  of  the  average  power  output  of  such  a  receiver  on  the 
mutual  intensity  functions  of  the  target-reflected  radiation  and  the 
local  oscillator  field  on  the  photodetector's  active  surface.  He  has 
found  that  the  dependence  is  through  only  the  real  part  of  the  product 
of  the  two  mutual  intensity  functions.  Since  the  mutual  intensity 
function  of  the  local  oscillator  is  real  and  constant  over  the  photo¬ 
detector  surface  for  normal-incidence,  uniform-intensity  plane-wave 
local  oscillator  beams,  it  follows  that  the  heterodyne  power  depends  on 
only  the  real  part  of  J0  in  such  cases  (assuming  the  photodetector  lies 
in  our  observation  plane ) . 

Figures  12  to  14  show  plots  of  the  real  part  of  JB(pj,Pi)  for  the 
cases  we  have  calculated.  Despite  the  apparently  discontinuous  char¬ 
acter  of  the  phase  of  the  mutual  intensity  as  a  function  of  d,  the  real 
part  shows  fairly  regular  behavior.  In  these  figures,  though,  it  is 
apparent  that  the  dependence  of  phase  on  0  becomes  much  stronger  as  R1 
is  increased.  The  dependence  of  phase  on  0  manifests  itself  as  the 
spreading  of  the  real  part  of  the  mutual  intensity  for  the  various 
values  of  0  plotted.  Figure  14  shows  the  maximum  spreading  observed. 
It  is  interesting  to  notice  that  in  figure  14  the  behavior  for  0  =  0 
closely  matches  the  behavior  for  0  =  tt;  similarly,  the  behavior  for  0  = 
n/4  matches  that  for  0  =  3tt/4.  It  thus  might  be  a  good  first  approxi¬ 
mation  to  characterize  the  real  part  of  the  mutual  intensity  as  follows. 
At  R1  =  0,  the  real  part  is  a  function  of  d  as  shown  in  figure  12  and 
has  circular  symmetry  about  Pj .  As  R1  increases,  the  real  part  remains 
similar  in  appearance,  but  the  symmetry  about  Pj  becomes  elliptical, 
with  the  major  axis  of  the  ellipse  of  symmetry  perpendicular  to  the 
radius  through  Pj .  The  eccentricity  of  the  ellipse  increases  with  R1 . 


7Dennis  (f.  McGuire,  Coherent  Detection  of  Partially  Coherent  Sources, 
Opt.  Lett,,  5,  2  (February  1980),  73-7 5. 
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Figure  14.  Mutual  intensity  versus  observation  point  separation 
(real  part) — fixed  point  3  cm  from  center. 


5.  DISCUSSION 


It  is  not  the  purpose  of  this  report  to  quantitatively  describe  the 
mutual  intensity  function,  but  rather  to  describe  how  this  function  can 
be  calculated.  We  have  demonstrated  that  these  calculations  can  be 
carried  out  within  the  limits  set  at  the  beginning  of  this  report. 
However,  as  was  stated  earlier,  these  calculations  are  only  a  first 
step.  The  final  objective  is  to  examine  the  utility  of  optical  heter¬ 
odyne  systems  in  laser  fuzing  systems.  To  this  end,  two  areas  of  con¬ 
tinued  research  are  outlined. 


The  first  area  concerns  a  problem  that  should  be  considered  in  any 
future  work,  which  is  that  a  laser  fuzing  system  generally  contains  a 
lens  system  in  the  receiver  to  collect  and  focus  the  received  field  onto 
a  photodetector.  The  mutual  intensity  propagation  equation  considered 
in  this  report  describes  free-space  propagation  of  light;  it  does  not 
cope  with  propagation  through  lenses.  Centered  collecting  optics  can  be 


dealt  with  by  calculating  the  mutual  intensity  at  the  entrance  pupil, 
mapping  the  entrance  pupil  response  to  the  exit  pupil  using  principles 
of  ray  optics,  and  then  propagating  the  mutual  intensity  at  the  exit 
pupil  onto  the  observation  plane  (photodetector).  The  last  step  is 
quite  problematic,  however,  because  the  mutual  intensity  at  the  exit 
pupil  (the  "target  plane"  for  the  second  propagation  calculation)  is  not 
necessarily  dependent  solely  on  |P.|-P2|.  At  present,  we  have  not  inves¬ 
tigated  how  this  should  be  handled. 

A  related  problem  occurs  in  considering  the  power  output  from  a 
photodetector  in  a  heterodyne  system.  In  this  case,  an  expression 
involving  the  mutual  intensity  at  the  observation  (detector)  plane  must 
be  integrated  over  the  detector  surface.  Again  the  function  to  be 
integrated  would  probably  not  be  a  function  of  |P}-P2I  alone.  The 
integrand,  however,  is  somewhat  simpler,  since  no  propagation,  and  hence 
no  path  length  differences,  are  involved. 

Both  of  the  above  calculations  would  be  simplified  considerably  if 
the  mutual  intensity  at  the  observation  plane  could  be  decomposed  into  a 
-urn  of  appropriate  radially  dependent  functions.  Lens  systems  could 
then  be  dealt  with  using  the  techniques  of  this  report.  (However,  edge 
effects,  which  we  have  not  considered,  may  be  too  prominent  in  this  case 
to  be  ignored. )  Also,  calculating  detector  power  output  would  be 
reduced  to  calculating  a  sum  of  integrals,  each  of  which  is  in  essence  a 
convolution  of  a  radially  dependent  function  with  the  detector  surface. 
The  first  proposed  area  of  future  research  is  to  investigate  how,  or  if, 
mutual  intensity  functions  calculated  using  the  techniques  of  this 
report  can  be  thus  decomposed. 

The  second  area  of  continued  research  is  to  develop  a  taxonomy  and  a 
catalog  of  responses  to  step-function  mutual  intensities.  This  would 


allow  a  parametric  characterization  of  the  mutual  intensity  at  an  obser¬ 
vation  plane  in  terms  of  the  cutoff  radius  rc  of  J  ( !  pi'  ~p2 1  » rc )  and 
geometrical  parameters.  Such  a  catalog  would  be  the  first  step  toward 
achieving  a  better  understanding  of  the  near-field  behavior  of  the 
mutual  intensity  function,  which  is  an  essential  feature  for  the  heter¬ 
odyne  system  analysis  because  of  the  relatively  short  detection  ranges 
involved  in  fuzing.  If  this  information  were  to  be  developed,  many 
expensive  intermediate  calculations  could  be  avoided. 
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APPENDIX  A. — FORTRAN  PROGRAM  TO  CALCULATE  PROPAGATION  OF  MUTUAL 
INTENSITY  FUNCTION 


APPENDIX  A 


C 

PROGRAM  RAYS 

C 

PURPOSE 

C 

THIS  PROGRAM 

CALCULATES  THE 

MUTUAL  INTENSITY  BETWEEN  TWO 

C 

POINTS,  ONE 

BEING  P1XED  AND 

THE  OTHER  MOVING  AWAY  ALONG 

c 

C' 

SELECTED  RAYS  AT  VARIOUS  ANGLES  PROM  THE  RADIUS. 

V. 

c 

USAGE 

c 

THIS  PROGRAM 

REQUIRES  A  SET 

OF  CONTROL  CARDS  ON  LOGICAL  UNIT 

c 

S  TO  DIRECT 

THE  CALCULATIONS.  THESE  CARDS  ARE  AS  FOLLOWS: 

c 

CARD  # 

CONTENTS 

FORMAT 

c 

I 

NPROB 

13 

c 

(2 

RANGE 

D1S.7 

c 

WAVE 

D1S.7 

c 

RSIGMA 

D15.7 

c 

3 

RCORR 

D1S.7 

c 

SIGMA2 

D15.  7 

c 

4 

NPOSN 

13 

c 

(5 

Kl 

D  IS.  7 

c 

BRR 

D1S.7 

c 

ICD 

14 

c 

6 

THSTEP 

D1S.7 

c 

THSTRT 

D1S.7 

c 

NTHSTP 

14 

c 

7 

DSTEP 

D1S.7 

c 

DSTART 

D1S.7 

c 

NDSTEP 

14  )  ) 

c 

IN  THE  ABOVE 

,  THE  SEQUENCE 

OF  CARDS  AFTER  CARD  1  IS  REPEATED 

c 

NPROB  TIMES 

(NPROB  IS  THE  NUMBER  OP  PROBLEMS).  AFTER  EACH 

c 

OCCURRENCE  OF  CARD  4  (ONCE 

IN  BACH  PROBLEM  SET),  CARDS  S 

c 

THROUGH  7  ARE  REPEATED  NPOSN  TIMES  ( NPOSN  IS  THE  NUMBER  OF 

c 

r 

SETS  OF  POSITIONS  FOR  OBSERVATION  POINTS). 

c 

THE  PROGRAM 

PERFORMS  THE  CALCULATION  OP  MUTUAL  INTENSITY  FOR 

c 

EACH  PROBLEM 

AND  EACH  POINT 

PAIR  AS  INDICATED  BY  THE  CONTROL 

c 

CARDS  ABOVE. 

THE  POSITION 

OF  THE  OBSERVATION  POINTS  ARE 

C  DETERMINED  BY  THE  CONTENTS  OF  CARDS  5-7  AS  FOLLOWS.  R1  IS  THE 

C  DISTANCE  OK  THE  FIXED  POINT  FROM  THE  CENTER  OF  THE  TARGET  SIGMA. 

C  IERR  AND  ICD  ARE  EXPLAINED  BELOW.)  THE  SECOND  POINT  IS  MOVED 

C  AWAY  FROM  THE  PIRST  ALONG  LINES  AT  VARIOUS  ANGLES  PROM  THE  RADIUS 

C  PASSING  THROUGH  THE  FIXED  POINT.  THERE  ARE  NTHSTP  OF  THESE 

C  LINES.  AT  ANGLES  THSTRT «  THSTRT ♦THSTEP,  THSTRT* 2*THSTEP,  ETC. 

C  ALONG  EACH  L I NE ,  THE  SECOND  POINT  IS  LOCATED  AT  VARIOUS  SEPA- 

C  RATIONS  FROM  THE  FIXED  POINT.  THERE  ARE  NDSTEP  POSITIONS  FOR 

C  THE  SECOND  POINT,  AT  SEPARATIONS  OF  DSTART ,  DSTART+DSTEP , 

C  DSTAKT+2*DSTEP,  ETC.  FROM  THE  FIXED  POINT. 

C 

C  THE  VARIABLES  ERR  AND  ICD  ARE,  RESPECTIVELY,  THE  REOUESTED 

C  MAXIMUM  RELATIVE  ERROR  OP  THE  NUMERICAL  CALCULATIONS  AND  A 

C  CODE  REQUESTING  THE  CALCULATIONS  DESIRED.  IP  ICD=1,  ONLY 

C  THE  RBSI DUAL  CONTRIBUTION  TO  THE  MUTUAL  INTENSITY  IS  CALCU- 

C  LATED,  OTHERWISE  THE  PULL  MUTUAL  INTENSITY  BETWEEN  EACH  PAIR 

C  OP  OBSERVATION  POINTS  IS  OBTAINED  (WITHOUT,  HOWEVER,  THE  EDGE 

C  EPFECTS  INCLUDED). 

C 

C  THE  REMAINING  VARIABLES  ABOVE  ARE  AS  FOLLOWS: 

C  RANGE  -  THE  OBSERVATION  DISTANCE  FROM  THE  TARGBT. 
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C  WAVE  -  THE  WAVBLBNCTH  OP  THB  ILLUMINATING  LIGHT. 

C  RSIOMA  -  THB  RADIUS  OP  THB  TARGET. 

C  RCORR  -  THE  PBRPBCT  CORRELATION  RADIUS. 

C  SIOMA2  -  THB  VARIANCE  OP  THE  SURPACB  HEIGHT  ROUGHNESS, 

C  SCALED  TO  THB  WAVELBNGTH. 

C  ALL  OP  THESB  QUANTITIES  ARB  INPUT  IN  MBTBRS,  EXCEPT  SIGMA2, 

C  WHICH  IS  DIMENSIONLESS. 

C 

C  NOTES 

C  THIS  PROGRAM  WAS  WRITTEN  TO  USE  A  VARIATION  OP  THE  ROUTINE 

C  I FBI X  OP  THE  HDL  PAGLOAD  LIBRARY  IN  WHICH  THB  NAME  POPX 

C  IS  ASSUMBD  POR  THE  SUBROUTINE  WHICH  COMPUTES  THE  INTEGRAND 

C  FOR  IPE1X.  THIS  WAS  DONE  POR  THE  SAKB  OP  RUNTIME  EFPICIBNCY. 

C 

C  IN  ESTIMATING  CPU  TIME  POR  THIS  PROGRAM  WHEN  THE  PULL  MUTUAL 

C  INTENSITY  IS  BEING  CALCULATED,  ALLOW  60  CPU  SECONDS  POP  EACH 

C  POINT  PAIR  WHEN  RUNNING  ON  THE  370/168  AFTER  BEING  COMPILED 

C  WITH  THE  FORTRAN  H-EXTBNDED  COMPILER  WITH  OPTIMIZATION  LEVEL  2. 

C 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMPLEX'*  16  JB  ,  JRES 
REAL* 8  JNMACfJRPN 

COMMON  /EVCN/PANGE, WAVER, 8RR/PROBLM/RS IGMA , RCORR, SIGMA2 

COMMON  /PARM/D,R1 , THETA 

DATA  Pl/3. 141S92653S89783238462643D0/ 

CALL  PGMASK 
C 

C  READ  NUMBER  OF  PROBLEMS  AND  LOOP  FOR  EACH  PROBLEM. 

C 

READ  (5,5000)  NPROB 
5000  FORMAT  (13) 

DO  30  VP=1 , NPROB 
C 

C  BEGIN  PROCESSING  -  READ  PROBLEM  PARAMETERS 
C 

READ  (5,5010)  RANGE, WAVE,RS1CMA, RCORR ,SIGMA2 
5010  FORMAT  (3D15.7) 

WAVE  IC=2  •  DO* PI  /WAVE 

WRITE  (6,6000)  RSIGMA , PANGS, RCORR, SIGMA2, WAVE 
6000  FORMAT  (•  TARGET  RADI US = • ,2PF8 .6 , »  CM,  TARGET  DISTANCE- • ,0PF5.0/ 

*  •  CORRELATION  RADI US=* , 2FF 1 0.8 , •  CM,  SURPACE  HEIGHT  VARIANCE—  * , 

*  0PF8. 4 , ' ,  WAVELENGTH  =  • , 6PF6 .2 , •  MICRONS*) 

C 

C  READ  NUMBER  UP  POSITION  PARAMETER  SETS  I  LOOP  FOR  EACH  ONE. 

C 

READ  (5,5000)  NPOSN 
DO  30  NPOS=l, NPOSN 

READ  (5,5020)  R1,ERR,I CD , THSTEP, THSTRT , NTHSTP ,DSTEP , DSTART , NDSTEP 
5020  FORMAT  (2D1S.7,I4) 

Kl=DABS ( R I ) 

C 

C  COMPUTE  ANGLE  OP  RAY  FROM  RADIUS  -  NORMALIZE  TO  BE  B8TWEBN  (0,PI). 

C 

DO  30  NTH=l  ,NTHSTP 

THETA=PI*(THSTBr*THSTBP*DPLOAT(NTH-l ) ) 

IF  ( THETA .GE .0. DO • AND .THETA. LE .PI )  GO  TO  20 
IF  ( THETA .LT.O. DO)  THETA«THETA*P I 
IF  (  THBTA.GT.PI )  THETA “THETA— PI 
GO  TO  10 
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C  STEP  OUT  ALONG  EACH  RAY  -  CALCULATE  GBOMBTRY 
C 

20  DO  30  ND  *1  , NDSTEP 

D=DSTART*DSTBP*DFLOAT( ND— 1 ) 

X2*Rl-D*DCOS(  THETA) 

Y2*D*DS1N (THETA) 

R2*DSORT( X2*X2*Y2*Y2) 

WRITE  (6,6010)  Rl , X2, Y2 ,D 

6010  FORMAT  ( • OOBSER VAT ION  POINTS  LOCATED  AT  (*,2PF6.3,*,  0.0  ),  (*, 

*  P6.3,  •  ,  •  ,66.3,  •  )  CM.;  SEPARATION  sSPi.a,'  CM.*) 

C 

C  CALL  GETANS  TO  OBTAIN  RESIDUAL  AND,  IF  DESIRED,  TOTAL  RESPONSE. 

C  IF  ICD-1,  ONLY  THB  RESIDUAL  IS  COMPUTBD,  AND  TUB  TOTAL  IS  SET 

C  TO  THIS  VALUES  OTHERWISE,  THE  TOTAL  RESPONSE  INCLUDES  THE  SHORT 

C  RANGE  CORRELATION  CONTRIBUTION. 

C 

CALL  GETANS ( Rl  , R2, JRBS ,JB,ICD, 1C HECK) 

C 

C  COMPUTE  MAGNITUDE  AND  PHASE  OF  TOTAL  RESPONSE  AND  PRINT  ANSWER. 

C 

JBMAG  =CDA  BS  (  J  B ) 

JBPH— 0.D0 

IF  (JBMAG.GT.O.DO)  JBPH-DATAN2< D IMAG ( JB ) ,DREAL( JB ) ) /PI 
WRITE  (6,6020)  JRBS,JB, JBMAG, JBPH 

6020  FORMAT  (•  JEES=( 1 , 1PD15 .8 , * , * ,D1S .8 , * ) ,  J B=( • ,DI5.8, • , » ,D1S .8 , 

*  •),  MAG—  * ,D15. 8, * ,  PHASE—* ,Dl 5.8 ) 

IF  ( ICHECX.LT. O)  WRITE  16,6030)  ICHECK 

6030  PORMAT  (SX,**«*  ICHECK  =• , 13,*  ***•) 

30  CONTINUE 
STOP 
END 

SUBROUTINE  GETAN S I R 1 , R 2 , JRBS , JB , I CD , I CH BCK ) 

C  PURPOSE 

C  THIS  SUBROUTINE  COMPUTES  THB  MUTUAL  INTENSITY  BETWBBN  TWO 

C  DATA  POINTS,  NEGLECTING  EDGE  EFFECTS.  IF  ICD=I »  ONLY  THE  RESIDUAL 

C  CONTRIBUTION  IS  CALCULATED  AND  JB  IS  SET  EQUAL  TO  JRES. 

C 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMPLEX* 16  EVJ, JB, JRBS,OSUBAC, INTGND 
EXTERNAL  INTGND 

COMMON  /EVCM/RANGB,WAVBK ,EKR/PROBLM/RS IGMA , RCORR , SIGMA2 
DATA  PI/3.141S92653589793238462643D0/ 

RES— DEXP ( -S IGMA 2 ) 

RMAX*Rl ♦ RS IGMA- RCORR 
SCALE*! WAVBXWRAMOB/ (2 . DO*PI ) ) **2 
C 

C  CALCULATE  SCALED  RESIDUAL  CORRELATION  CONTRIBUTION 
C 

JRESsRBS* EVJ ( Rl , RS IGMA) *DCON JC ( BV J ( R2 , RSIOMA )  ) 

JB*JRES 

C 

C  IP  REQUESTED,  1NCLUDB  THE  CONTRIBUTION  PROM  SHORT  RANGE 
C  CORRELATION,  EXCLUDING  EDGE  EFFECTS. 

C 

IP  (ICD.NE.i)  JB=JRBS*< 1 .DO-RES ) *OSUBAC (0. DO, RNAX, ERR ,NPTS, 

*  ICHECK, RELERR, INTGND) 

J  B*SC  ALE  *  JB 

J  RES* SCALE* JRES 

RETURN 

END 
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COMPLEX  FUNCTION  I NTGND* 16 1  RAD IUS ) 

C  PURPOSE 

C  THIS  FUNCTION  EVALUATES  THE  INTEGRAND  OP  THE  OUTERMOST  INTEGRAL 

C  FOR  COMPUTING  THE  SHORT  RANGE  CORRELATION  CONTRIBUTION  TO  THB 

C  MUTUAL  INTENSITY  FUNCTION. 

C 

IMPLICIT  HEAL*8  JA-H,0-Z> 

COMPLEX* 16  EVJ ,  EVTH 

COMMON  /PROBLM/RSIGMA , RCORR/EVCM/RANGE t  WAVER 

COMMON  /PARM/D, R1 , THBTAR 

DATA  PI/3 .14 15926S3589793238462643D0/ 

C 

C  EVALUATE  FIRST  FACTOR  OF  INTBGRAND  -  RESIDUAL  PROM  A 
C  SMALL  CIRCLE. 

C 

INTGND*RADIUS*BVJ I  RADIUS, RCORR) 

C 

C  EVALUATE  SECOND  FACTOR  OF  INTBGRAND  -  AN  INTEGRAL  IN  THETA 
C 

IF  ( RADIUS. GT. RSIGMA— RCORE—R 1 )  GO  TO  10 

C 

C  TBST  FOR  DEGENBRATE  CASB 
C 

THBTAl=PI 

IF  (D.EO.O.DO)  GO  TO  SO 
C 

C  INTEGRAL  IN  THETA  IS  OVBR  A  PULL  CIRCLE  -  INTEGRAND  IS  TWICE 

C  THE  INTEGRAL  FROM  0  TO  PI  IN  THETA,  TIMES  THE  PIRST  FACTOR. 

C 

1NTGND=I NTGND *2 .D04BVTH ( RADIUS , 1 .DO , -I . DO ) 

RETURN 

C 

C  INTEGRAL  IN  THETA  IS  OVBR  LESS  THAN  A  FULL  CIRCLE  -  SELBCT 

C  THB  RIGHT  BHEAXUP  OP  THIS  INTBGRAL  AND  COMPUTE  INTBGRAND 

C  AS  THE  PRODUCT  OP  THIS  INTBGRAL  WITH  THE  FIRST  FACTOR. 

C  PIRST,  FIND  THB  ENDPOINTS  OP  THB  INTEGRATION  IN  THETA.  THB 
C  ENDPOINTS  ARE  THE TAR— THETA 1  AND  THETAR+THBTA1 • 

C 

10  X1*(RADIUS*RADIUS4R1*R1-(RSIGMA-RCORR)**2)/(2.DO*R1*RADIUS) 
THETA1 =DARCOS (X 1 I 
IP  (D.EQ.O.DO)  GO  TO  50 
II^COSl  THETAR I 

XA*X1*XR-DSQRT( ( 1 .D0-X1 *X1 >• ( 1 .DO— XR*XR ) 1 
XB*X1«XR4DSQRT( (1 .DO— XI *X 1)* ( 1 .D0-XR4XR >) 

IP  | THETAR— THETA1 .LT.O .DO .AND.THET AR4TH BTAl .LE« PI 1  GO  TO  20 
IF  ITHETAR-THETA1.GE.0. D0.AND.THETAR4THBTA1.GT.PI)  GO  TO  30 
IF  (THBTAR-THETA1.LT.0. DO. AND.THETAR4THETA1.GT.PI )  GO  TO  40 
C 

C  CASE  1  -  0  <=  THBTAR— THETA 1  <*  THBTAR 4THETA1  <*  PI 
C 

1NTGND=I NTGND4BVTH ( RAD I US ,XB ,XA) 

RETURN 

C 

C  CASE  2  -  TBBTAE-THBTA1  <  0  <=  THBTAR4THBTA1  <*  PI 
C 

20  INTGND*INTOND*(2.D0*BVTB(RADIUS,l.D0,XB)4BVTH(RADIUS,XB,XA) ) 
RETURN 

c 

C  CASE  3  -  0  <*  THBTAR— THETA 1  <*  PI  <  THBTAR4THETA 1 
C 
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30  INTGND=I NTGND*  ( EVTH ( RADI U S  ,XB ,XA ) +2  .D0*BVTH (RADI US, XA ,— 1 »D0 ) ) 

RETURN 
C 

C  CASE  4  -  THETAR— THETAl  <  0  <  PI  <  THETAR ♦THETAl 
C 

40  I NTGND*I NTGND* ( 2. DO*BVTH ( BAD I US , 1 .DO , XB  )♦ 

*  EVTH ( RAD IUS, XB ,XA ) 42  .DO* EVTH ( BAD I US , XA,  — 1 .DO ) ) 

RETURN 

C 

C  DBGBNERATE  CASE  -  INTBCRAL  IN  THETA  CAN  BE  COMPUTED  DIRECTLY 
C 

50  1NTGND=INTGND*2.D0*THETA1*CDEXP(DCMPLX( 0. DO, -WAVE!* 

*  DSORT( RANGE**2*RADIUS**2) ) )/( RANGE**2 ♦RADIUS** 2 ) 

RETURN 

END 

COMPLEX  FUNCTION  BVJ*16 (R, RADIUS) 

C  PURPOSE 

C  THIS  SUBBOUT INB  CALCULATES  THE  MUTUAL  INTENSITY  FUNCTION  FOR 

C  POINTS  SEPARATED  BY  A  DISTANCE  »R •  FROM  A  TARGET  OF  RADIUS 

C  'RADIUS'. 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX*  1 6  D,  Z 

COMMON  /EVCM/ RANGE ,WAVEK , ERR/ EVJ CM/RAD , OFFSET 
COMMON  /CONTRL/ISWTCH 
LOGICAL  LONG 

DATA  PI/3.14 15926535 89 793238462643D0/ 

C 

P(X) =( I .DO— 2. DO* ( 1 . DO— I  .2Dl*(t .D0-3 .D1 / ( X*X ) ) /( X*X>  >/<X*X) )/X 
O 1X1  =  11. DO-6. DO* ( 1. DO-2. Dl*( 1 .DO-4  .2D1 / (X*X ) ) /( X*X) ) / ( X*X ) ) / (X*X ) 

C 

RAD=RAD1US 

OFFSET*R 

Z=CDEXP(DCMPLX(O.DO,WAVEK*RANGE> I 
LONG RAD IUS— R ) .LT • 1 • D— 2* RANGE 
SXD=(  (  RAD  IUS— R)  /RANGB)  **2 
A=WAVEX*RANGE*(DSORT( 1 .D0*SKD)-1 .DO) 

IF  (LONG)  A=VAVBK* RANGE *SKD*( l.DO-SKD*( 1. DO-SKD/2. DO ) /4 .DO) /2.D0 
SXD»WAVBK*RANGE*DSORT ( 1 • DO+SKD) 

c 

C  CALCULATE  CONTRIBUTION  FROM  INSIDE  TARGBT  CYLINDER  (IF  ANY). 

C 

D~(0.D0,0.D0) 

IF  ( R.LT. RADIUS )  D*2.D0 *PI*Z*DCMPLX ( 

*  G( WAVEK*RANGE) ♦F(SKD)*DSIN(A)-G(SKD) *DCOS(  A)  , 

*  P ( WA VEX*RANGB )— F (SXD )«DCOS (A )— G ( SKD ) *DSIN ( A )  ) 

C 

C  CALCULATE  RBMAINDBR  OF  CONTRIBUTIONS 
C 

B  =  (( RADI US*R) /RANGE)** 2 

IF  (LONG)  Bx*AVEK*RANGE*B*(l.DO-B*(l.DO— B/2.D0)/4.D0)/2.D0 
IF  (.NOT. LONG)  BsWAVBK* RANGE* ( DSORT ( 1 .DO*B) — 1 .DO) 

ISWTCH=l 

IBR*1 

CALL  IFBIXiEREAL, BINAG, A, B, 100000 ,  . 1 D0*BRR ,0 , 0 , 1ER ) 

IF  (IBR.E0.2)  GO  TO  10 

EVJ =D*2 . DO* Z* DCNPLX ( ERE AL ,S I MAG ) 

RETURN 

10  WRITE  (6,6010)  R, RADIUS, A 

6010  FORMAT  (•  IFEIX  DID  NOT  CONVERGB  -  CALLED  FROM  EVJ.'/ 
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*  •  VARIABLES  ARB :  R= • , 1 PDI 3.6, • ,  RAD I  US* • * D 13 .6 * 

*  »,  A=',D13.6/»  PROGRAM  STOPPED' I 
STOP 

END 

COMPLEX  FUNCTION  BVTH* 16( OFFSET, COSO ,COSI ) 

C  PURPOSE 

C  THIS  SUBROUTINE  COMPUTES  THB  MUTUAL  INTENSITY  DUB  TO  SHORT 

C  RANGE  CORRELATION  PROM  A  TARGET*  WITH  ONE  OBSERVATION  POINT 

C  SEPARATED  BY  'D<  FROM  THE  OTHER  ONE. 

C 

IMPLICIT  KEAL*8  (A-H.O-2) 

LOGICAL  LONG 

COMMON  /EVCM/RANGB*WAVEX*BRR/EVTHCM/A*B,BLESSA 
COMMON  /CONTRL/1SWTCH/PARM/D 
LONG-( (O-OFFSET) /RANGE t **2.LT.I.D-3 
As*AVEK*DSQRT(RANGB**2+(D-OFPSET  >**2) 
B3WAVEK*DSORT(RANGE**2+(D+OFPSET)**2> 

BLESS A^B-A 
IF  (LONG)  BLESSAs 

*  2.D0*WAVEK*D*OFFSBT*( l.D0-(OPPSET**2«D**2)* 

*  ( I .DO  — ( OFFSET** 2 AD** 2 ) /RANGE** 2 ) / ( 2 .DO  *RANOE  **2 )  — 

*  (OFPSBT**2-D**2)**2/(8.DO*RANGE**4))/RANGB 
C 

C  TBST  IF  NEAR  TOP  OP  INTERVAL  ( A  *  B ) 

C 

IP  (COSO.LT.O.DO)  GO  TO  30 
C 

C  NEAR  BOTTOM  END  -  TEST  IP  LONG  RANGE  APPROXIMATION  SHOULD  BE  USBD 
C 

IF  (LONG )  GO  TO  10 
C 

C  USB  SHORT  RANGE  EQUATIONS* 

C 

X0«WAV8X*D60RT( RANGB**2*D**2*OFFSET**2-2. D0*D*OFPSET*COS0 )-A 
X1«WAVEK*DSQRT( RANGE* *2 ♦D**2*QFFSBT**2— 2. DO*D*OPFSBT*COS1 )-A 
GO  TO  20 
C 

C  USB  LONG  RANGE  BOUATIONS  WITH  INTERVAL  NOT  NEAR  TOP  END  OF  (A»B) 

C 

1 0  X0=*D*OPPSBT* ( 1 . DO-COSO ) 

X l«D*OFFSET* ( 1 • DO-COS  I ) 

X*( D—OPFSBT) **2 

X0*WATBX*X0* ( 1 .DO— ( X+XO ) • ( 1.D0— ( X*X0 )/R ANOE**2 ) /( 2 .DO*  RANGE**  2 )• 

*  X* ( X+2.D0*X0 ) / ( 8 .DO* RANGE* *4 ) ) /RANGE 

X1*VAVBK*X1* ( 1 .DO— ( X+Xt )• ( l.D0-( X*X 1 > /RANGE* *2 ) /( 2. DO* RANGE** 2)  • 

*  X*(X*2.D0*X1)/(8.D0*RANCB**4) )/ RANGE 

20  X0«DMIN1 (BLBSSA*DMAXi 1X0*0.00)) 

X 1*DMAX1 ( 0. DO. DM INI ( XI  * BLESSA ) ) 

ISNTCEx2 

IBR'l 

CALL  IFB IX( ANSR *ANS1* XO * XI* 100000* . 1D0*BRR*0*0* 1BR) 

IF  ( IER.EQ. 2)  GO  TO  60 

BVTH«2.D0**AVEX**2*DCMPLX( ANSR.-AMS1) *CDEXP (DCMPLX ( 0 .DO , -A ) ) 
RETURN 
C 

C  NEAR  TOP  END  -  TBST  IF  LONG  RANGE  APPROXIMATION  SHOULD  BB  USBD. 

C 

30  IF  (LONG)  GO  TO  40 

XO*RAVBK*DSORT(RANGB**2*D**2*OPFSBT**2-2.DO*D*OFFSBT*COSO)-B 
Xt*VAVBK*DSORT( RANGB**2*D**2*OFFSBT**2-2. DO*D*OFPSET*COS1 ) -B 
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GO  TO  50 
C 

C  USB  LONG  RANGB  EQUATIONS  WITH  INTBKVAL  NBAS  TOP  END  OP  (A,B) 

C 

40  X0=—D*OPPSET*  1 1  «DO-*COSO) 

Xl=-D*OPFSET*(I .DO *008 I I 
X= I D ♦OFF SET 1**2 

X0-WAVBK*X0*(1 .D0-U*X0>* ( 1 .D0-( X*XO> /BAMOB**2 1 / 1  2 .D0*RANGB**2 >- 

*  X* (X*2 .D0*X0)/( 8. DO* RANGE* *4) ) /RANGB 
Xls*AVEX*Xl*(l.D0-|X4Xl I* (l.DO-(X*Xl I /RANOB**2 ) / t 2 .D0*BANGB**2 >- 

*  X*(X*2.1>0*xn/(8.D0*RANGB**4)  ) /RANGE 

SO  X03DMAX1(-BLESSA,DMIN1(X0,0.D0|) 

Xl*DMINl (O.DOyDMAXl (XI ,-BLESSA) I 
ISWTCH=3 
I  ER=1 

CALL  1PBIX(ANSB,ANSI,X0,X1,  100000,.  IDO*  ERR,  0,0,  IER) 

IP  (IER.EQ.2)  CO  TO  60 

EVTH*2.D0*WAVLX*«2*DCMPLX (ANSR,-AM8I ) *CDCXP( DCMPLXI O.DO ,*•  > I 
RETURN 

60  PRINT  *, OFFSET, COS0,COSi, A, B,BLBSSA, ISWTCH, XO, II 

STOP 
END 

SUBROUTINE  FOPX|XtY) 

IMPLICIT  RBAL*M  (A-H,0-X) 

COMMON  /EVCM/RANGE tWAVEK/BVJCM/R |S 
COMMON  /EVTHCM/ A, B, BL8SSA /CON TRL/ ISWTCH 
COMMON  /CMSPLT/AO 
IOPT— 0 

GOTO  ( 10  ,20,30) .ISWTCH 
STOP  1 
C 

C  INTEGRAND  FOR  EVJ 
C 

10  Y*=X/WAVBI 

Y=(Y*Y*2.DU*RANGE*Y*S*S-R*H|/(2.D0*S*DSQBT(Y*Y*2.D0*RANGB*TH 
IE  ( DABS (Y ) .GT. 1 . DO )  Y=DS IGN { 1 .DO , Y J 
Y=DARCOS ( Y ) / ( X*WAV EX* RANGE  I 
GO  TO  40 
C 

C  INTEGRAND  POR  EVTH 
C 

20  Y=l.DO/( (X*A)*DSORT(X*(BLESSA-X)*(X*A*B)*(X*2.D0*A) ) ) 

GO  TO  40 
C 

C  INTEGRAND  FOR  EVTH  -  TOP  OF  (A,B) 

C 

30  Y^l.DO/l ( X*B)*DSQRT(-X*(X+BLESSA>* (X+AEB) *(X*2.D0*B) ) ) 

40  IF  (IOPT.E0.il  Y=Y*( 1 .D 0*DCOS( X-AO )) 

IF  (IOPT.EO.2!  Y=Y*( 1 . DO+DSI N ( X-AO ) ) 

IP  (IOPT.EQ.3)  Y=Y*DCOS ( X  I 
IP  ( IOPT.EO.4)  Y=Y*DSIN(X) 

RETURN 

BNTRY  FCOSXA(X  »  Y  ) 

I  OPT* 1 

GO  TO  ( 10,20,30) , ISWTCH 
STOP  1 

ENTRY  PS INXA (X , Y I 
IOPT=2 

GO  TO  I  10, 20, 30 }, ISWTCH 
STOP  1 


63 


DISTRIBUTION 


ADMINISTRATOR 

DEFENSE  TECHNICAL  INFORMATION  CENTER 
ATTN  DTIC-DDA  (12  COPIES) 

CAMERON  STATION,  BUILDING  5 
ALEXANDRIA,  VA  22314 

COMMANDER 

US  ARMY  RSCH  S  STD  GP  (EUR) 

ATTN  CHIEF,  PHYSICS  £  MATH  BRANCH 
FPO  NEW  YORK  09510 

COMMANDER 

US  ARMY  MISSILE  S  MUNITIONS 
CENTER  £  SCHOOL 
ATTN  ATSK-CTD-F 
REDSTONE  ARSENAL,  AL  3 5809 

DIRECTOR 

US  ARMY  MATERIEL  SYSTEMS  ANALYSIS  ACTIVITY 
ATTN  DRXSY-MP 

ABERDEEN  PROVING  GROUND,  MD  21005 
DIRECTOR 

US  ARMY  BALLISTIC  RESEARCH  LABORATORY 
ATTN  DRDAR-TSB-S  (STINFO) 

ABERDEEN  PROVING  GROUND,  MD  21005 

TELEDYNE  BROWN  ENGINEERING 
CUMMINGS  RESEARCH  PARK 
ATTN  DR.  MELVIN  L.  PRICE,  MS-44 
HUNTSVILLE,  AL  35807 

ENGINEERING  SOCIETIES  LIBRARY 
ATTN  ACQUISITIONS  DEPARTMENT 
345  E.  47TH  STREET 
NEW  YORK,  NY  10017 

HQ,  USAF/SAMI 
WASHINGTON,  DC  20330 

US  ARMY  ELECTRONICS  TECHNOLOGY 
fi  DEVICES  LABORATORY 
ATTN  DELET-DD 
FT  MONMOUTH,  NJ  07703 

COMMANDER 

US  ARMY  ARMAMENT  RESEARCH 
fi  DEVELOPMENT  COMMAND 
ATTN  DROAR-SE,  SYSTEMS  EVALUATION 
OFFICE,  LTC  GRADY  COOK 
ATTN  DRDAR-LCN-C ,  G.  TAYLOR 
ATTN  DRDAR-AS,  R.  ANDREJZOVICS 
DOVER,  NJ  07801 


DIRECTOR 

NIGHT  VISION  S  ELECTRO-OPTICS 
LABORATORY 

ATTN  DELNV-V  (2  COPIES) 

ATTN  DELNV-R,  R.  BUSER 
FT  BELVOIR,  VA  22060 

COMMANDER/DIRECTOR 

ATMOSPHERIC  SCIENCES  LABORATORY 

ATTN  DELAS-EO,  F.  NILES 

WHITE  SANDS  MISSILE  RANGE,  NM  88002 

COMMANDER 

US  ARMY  MISSLE  COMMAND 
ATTN  DRSMI-REM,  H.  ANDERSON 
REDSTONE  ARSENAL,  AL  35898 

ENVIRONMENTAL  RESEARCH  INSTITUTE 
OF  MICHIGAN 
PO  BOX  8618 
ATTN  IRIA  LIBRARY 
ANN  ARBOR,  MI  48107 

COMMANDER 

NAVAL  WEAPONS  CENTER 

ATTN  CODE  3311,  K.  BULLOCK  (4  COPIES) 
ATTN  CODE  3311,  J.  CRISLER 
CHINA  LAKE,  CA  93555 

COMMANDER 

AFATL 

ATTN  L.  STABLES 

EGLIN  AIR  FORCE  BASE,  FL  32542 

COMMANDER 

AFAL/TEO 

ATTN  R.  HARRIS 

WRIGHT-P.*TTERSON  AF  BASE,  OH  45433 
MOTOROLA,  INC 

GOVERNMENT  ELECTRONICS  DIVISON 
ATTN  A.  GARAS  (2  COPIES) 

PO  BOX  1417 

8201  E.  MCDOWEL  RD 

SCOTTSDALE,  AZ  85252 

SANTA  BARBARA  RESEARCH  CENTER 
ATTN  N.  RIGBY  (2  COPIES) 

75  COROMAR  DRIVE 
GOLETA,  CA  93117 

NATIONAL  BUREAU  OF  STANDARDS 
ATTN  THEODORE  H.  HOPP  (10  COPIES) 
ROUTE  270  fi  QUINCE  ORCHARD  RD 
GAITHERSBURG,  MD  20760 


n 


Fi 


Fi 


rl 


» 


» 


9 


9 


65 


DISTRIBUTION  (Cont'd) 


US  ARMY  ELECTRONICS  RESEARCH 
&  DEVELOPMENT  COMMAND 
ATTN  TECHNICAL  DIRECTOR,  DRDEL-CT 

HARRY  DIAMOND  LABORATORIES 

ATTN  CO/TD/TSO/DI VISION  DIRECTORS 

ATTN  RECORD  COPY,  81200 

ATTN  HDL  LIBRARY,  81100  (2  COPIES) 

ATTN  HDL  LIBRARY,  81100  (WOODBRIDGE) 

ATTN  TECHNICAL  REPORTS  BRANCH,  81300 

ATTN  CHAIRMAN,  EDITORIAL  COMMITTEE 

ATTN  LEGAL  OFFICE,  97000 

ATTN  MORRISON,  R. ,  13500  ( GIDEP ) 

ATTN  CHIEF,  13000 

ATTN  CHIEF,  11000 

ATTN  CHIEF,  13300 

ATTN  CHIEF,  11400 

ATTN  COX,  L. ,  00211 

ATTN  HUMPHREY,  R. ,  13300 

ATTN  GRIFFIN,  J.  R. ,  21100 

ATTN  HATTER Y,  W.  V.,  13300 

ATTN  SANN,  K.  H. ,  15000 

ATTN  PEPERONE,  S. ,  36000 

ATTN  DOBRIANSKY,  B. ,  13500  (2  COPIES) 

ATTN  LAN HAM,  C. ,  00213 

ATTN  GIGLIO,  D. ,  15300 

ATTN  GEIPE,  T. ,  22100 

ATTN  TUTTLE,  J. ,  21400 

ATTN  VANDERWALL,  J. ,  22800 

ATTN  CONNER,  M.,  15200 

ATTN  TOZZI,  L. ,  34000 

ATTN  SOLN,  J.  R. ,  22300 

ATTN  BRANDT,  H. ,  22300 

ATTN  MCGUIRE,  D. ,  13300  (10  COPIES) 
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